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Final  Report:  PLASTIC  DEFORMATION  OF  GRANULAR  MATERIALS 

E.  Bruce  Pitman* 

April  1,  1988  -  March  31,  1989 

INTRODUCTION 

The  deformation  of  granular  materials  occurs  in  a  variety  of  applications:  soil  dynamics; 
avalanch  flo\vs:  grain  flow  in  bins.  Our  goal  in  this  project  is  to  develop  a  more  complete 
understanding  of  the  mechanics  of  granular  flow.  The  project  consists  of  two  distinct 
but  related  phases.  One  phase  is  a  study  of  the  mathematical  structure  of  constitutive 
relations  used  to  model  the  granular  medium.  The  research  here  is  primarily  concerned 
with  the  stability  and  well  posedness  of  the  evolution  equations  governing  flow.  Most  of  our 
efforts  have  concentrated  on  the  Critical  State  Theory  of  Soil  Mechanics,  a  mathematically 
attractive  theor\-  which  has  some  success  in  modeling  soil  deformations.  In  particular,  the 
Critical  State  Theory  is  reasonably  successful  at  modelling  the  deformation  of  so  called 
■‘Cam  clay",  a  type  of  clay  tested  by  the  Cambridge  soil  mechanics  group  in  the  late 
l950's.  Critical  State  Theory  applied  to  “Cam  clay”  is,  therefore,  a  useful  model  for 
soil  dynamics.  Complementing  this  theoretical  study  is  a  computationed  phase,  consisting 
of  numerical  simulation  of  cannonical  granular  flows.  In  Section  1  we  briefly  review  the 
equations  governing  granular  flow;  this  section  includes  very  preliminary-  ideas  regarding  a 
fully  non-linear  constitutive  relation.  In  Section  2  we  describe  the  results  on  the  stability 
of  flow.  In  Section  3  we  describe  work  on  the  computational  phase  of  tliis  project.  A 
review  of  many  of  these  issues  appeeirs  in  [8]. 

§1.  EQUATIONS  OF  MOTION 

(a)  Critical  State  Theory.  The  fundamental  equations  governing  the  deformation  of  a 
continuum  are  the  balance  laws  of  mass  and  momentum  : 

( 1.1 )  d,p  -t-  pV  ■  V  =  0 

(1.2)  ;'d.v  +  Vr  =  6 

Here  p  is  the  density  of  the  material  at  the  point  r  =  (xj,  xa,  ^3 )  and  time  t,  v  =  (ly .  V2.  V3  ) 
is  the  velocity,  and  T  is  the  syn.n.'  ti stress  tensor.  Body  forces  are  represented  by  h. 
cind  di  =  -I-  i’  •  V  is  the  convec-iv.-  derivitive.  We  employ  the  summation  convention 

throughout  this  report. 

Equations  1.1  -  1.2  represent  f  ,  i.a'ions  for  the  ten  unknowns  p.i\.T,y  In  ord.r 
close  the  system,  we  need  to  po--  ■  ■  •.'titutive  laws  relating  the  density  and  velocities 

‘Institute  for  .Mathematics  and  its  .-Xjii.li'  .i'.  .::'  I  :u\ersity  of  Minnesota.  Minneapolis,  Minnesota  .oS-tSS 
and  Department  of  Mathematics.  Nev^  .I>  rs4  \  Invctiite  of  Technology.  Newark.  New  Jersey  07102  This 
research  has  been  supported  by  the  .Air  I<irr»-  (iffire  of  Scientific  Research  under  grant  8!<-0182 
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to  the  stresses.  In  [9,11],  Pitman  and  Schaeffer  employed  the  theory  of  Critical  State  Soil 
Mechemics  [3,12]  to  provide  these  relations.  If  we  regard  compressive  stresses  as  positive 
and  denote  the  streiin-rate  tensor  as  Vij  =  —j(djV,  +  d,Vj),  these  constitutive  laws  may  be 


written: 

(1.3) 

0(  T,  p)  =  0 

(1.4) 

►,  ■  •• 

II 

Equation  1.3  is  the  Plastic  Yield  Condition  and  Equation  1.4  is  the  Associated  Flow  Rule. 
In  particular.  Equation  1.4  states  that  the  normal  to  the  Yield  Surfaice  (p  =  0  is  proportional 
to  the  strain-rate  tensor.  Note,  however,  that  the  constant  of  proportionality,  /i,  is  not 
specified  a  priori  (as  in  elasticity)  but  must  be  determined  as  part  of  the  solution  procedure. 
Different  versions  of  the  theory  may  be  examined  by  altering  the  function  d>.  It  is  usually 
required  that  <5  be  a  convex  function  of  T  and  monotone  in  p.  This  monotonicity  generates 
a  nested  sequence  of  yield  surfaces  in  stress  space;  such  a  picture  lies  at  the  heart  of  the 
Critical  State  Theory.  (Remark:  For  some  materials  like  dry  sand.  0  may  be  taken  to  have 
the  form 

0(7,  p)  =  >p(tr(T),  jdei;(r)|)  -  p^ . 

Here.  tr{T)  denotes  the  trace  of  T;  |der(T)i  is  the  Euclidean  norm  of  the  deviator  of  T. 
where 

*1.(1)  = 

and  is  a  sm.all  parameter  which  measures  compressibility,  t}-pically  (3  10“^  —  lO"^.) 

The  Associated  Flow  Rule  1.4  is  not  universally  accepted  in  plasticity  research,  espe¬ 
cially  met2d  plasticity  [2];  associative  flow  rules  were  introduced  in  plasticity  theory  by 
analogy  with  work  in  elasticity.  See  subsection  (b)  for  a  different  approach.  Standau-d 
non- Associative  theories  postulate  a  “flow  potential”  d'  and  the  flow  rule  We 

consider  only  the  Associative  theory  here. 

We  claim  that  the  system  of  equations  1.1  -  1.4  is  closed.  There  are  now  11  unknowns. 
p  being  added  to  the  previous  list;  there  are  also  11  equations  with  which  to  determine 
these  variables.  However,  Equations  1.3-1. 4  are  not  evolution  equations,  and  the  mixed 
JifTerential/aJgebraic  character  the  system  is  difficult  to  analyze.  Considerable  sim¬ 

plification  results  if  we  solve  1.3  1  «  for  T  and  p  as  functions  of  1'  and  p  and  write 

(1.5)  T  =  T(V,p). 

We  remark  in  Equa'icn  ?  F.  •  ;.•  ..  i^n^jutieci  as  a  function  of  1’.  is  homogeneous 
of  degree  zero;  this  homogeneity  i-  •  .u  teristic  of  a  rate  independent  theory. 

Before  finishing  this  section,  v..-  .-ome  terminology-  which  shall  be  used  subse¬ 

quently.  Let  e,  denote  the  eigenv;;i  ;ov  ..f  the  strain-rate  tensor  1’.  By  fully  three  dimen¬ 
sional  flow  we  mean  flows  for  wliirli  f,  ^  0  for  i  =  1.2.3.  Two  dimensional  flow.  then, 
meeins  either  flows  for  which  one  of  the  c,  =  0  or  flow  in  only  two  space  dimensions  xj.  j-t. 
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(b)  Alternative  Constitutive  Relations. 

The  Critical  State  theory  predicts  the  initial  deformation  of  clays  and  soils  reasonably 
well.  It  has  difficulties,  however,  with  cyclic  loading  paths  and  with  sandy  materials.  In 
addition,  when  the  Critical  State  theory-  predicts  ill  posedness  there  is  no  mecheinism  in 
the  theory  which  allows  for  recovery.  (We  skirt  this  difficulty  by  adding  kinetic  theory 
terms;  see  Section  3).  This  defect  is  shared  by  most  standard  deformation  theory  models. 
Prehminary  work  has  begun  to  investigate  a  class  of  “Yield- Vert  ex”  theories  of  the  type 
introduced  by  ChristofFersen  and  Hutchinson  [2]  for  polycrystalline  metals. 

The  fundamental  relation  in  the  Yield- Vertex  models  is  a  fully  non-linear  stress  strain- 
rate  relation  of  the  form: 

(1.6)  V',  = 

Here.  T  denotes  the  Jauman,  or  co-rotating,  time  derivitive  of  the  stress  T.  For  a  rate- 
independent  theory,  "k  must  be  homogeneous  of  degree  1  in  T;  symmetiy-  conditions  place 
other  constraints  on  the  form  of  (,1-6).  Modelling  is  required  in  order  to  correctly  account 
for  proportional  aind  non-proportional  loading  and  for  unloading  in  granular  materials. 

To  summarize  the  work  to  date.  Schcieffer  emd  Shecirer  have  begun  an  analysis  of  the 
structure  of  the  Yield- Vertex  theory.  In  state  space,  the  non-h\-perbolic  region  appears  to 
be  of  finite  extent;  it  is  possible  for  the  linearized  equations  to  become  ill  posed,  but  for 
the  system  to  remain  non-linearly  stable.  Pitman  has  performed  munerical  tests  on  model 
problems  in  order  to  gain  an  understanding  of  the  dynamics  predicted  by  the  theory.  A 
full  investigation  of  the  theory  is  planned. 

§2.  INSTABILITY 

We  may  substitute  Equation  1.5  into  1.1-1. 2  to  obtain  a  system  of  four  time  dependent 
partial  differential  equations  in  the  variables  p, Ui,f2,i’3  given  by: 

(2.1)  d,p  +  pd,Vj  =  0 

(2.2)  pd,i-.  +  (^)d,P  -  =  »• 

Here  we  have  dropped  body  forces  from  2  2.  Equations  2. 1-2.2  formally  resemble  the 
Navier-Stokes  equations  for  a  viscous  compressible  fluid.  However,  this  apparent  similarity 
is  illusary.  The  primaiy-  difference  Ixuween  the  fluid  equations  and  the  system  above  is  in 
the  nature  of  the  viscous  dissipation  term  In  the  Navier-Stokes  system,  the  symbol 
matrix  associated  to  this  dissipation  is  negative-definite;  in  Equations  2. 1-2.2,  the  symbol 
is  only  semi-definite.  The  non-d'  fiuiteness  of  this  symbol  is  the  origin  of  many  of 
^ornnlex  s-;acrures  in  granulm  fiwA-. 

Schaeffer  and  Pitman  [9,11]  exaininffi  the  linear  well  posedness  of  Equations  2. 1-2.2. 
In  particu'ar.  they  found  that  fully  three  dimensional  flows  are  well  posed.  However  two 
dimensional  flows  may  be  ill  joo^ed  with  a  growth  rate  of  Ollf').  where  f  is  the  Fourier 
tariable  dual  to  r.  Our  studies  in  this  phase  of  the  project  have  examined  the  stability  of 
the  Critical  State  equations  [10,13].  In  this  section  we  summarize  those  results. 
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For  many  typical  granular  flows,  there  exist  two  well  separated  time  scales  [9].  It  is 
possible  to  exploit  these  separate  scales  and  incorporate  into  our  analysis  the  effects  of 
terms  not  included  in  [9,11].  To  this  end,  define  a  new  time  f  =  cf  and  a  similarly  scaled 
velocity  v  =  Substitute  into  Equations  2. 1-2.2  and  recall  the  homogeneity  of  T{\\p) 

to  find  (dropping  the  tilde) 

(2.3)  dtp  +  dj{pvj)  =  d 


(2.4)  i'^pdtv,  + 

Now  make  the  quasi-dynamic  approximation  by  setting  e  =  0.  Linearize  the  equations 
about  a  steady,  homogeneous  solution  poi^’o-  In  the  linearization,  four  terms  arise  from 
the  derivitive  djipvj);  by  appropriate  change  of  coordinate  systems  we  can  drop  all  of 
these  terms  except  podjVj.  See  [13]  for  details.  We  derive,  then,  a  system  of  one  evolution 
equation  for  p  coupled  to  three  apparently  elliptic  equations  for  the  vp. 

{2.3*)  dtP  +  PodjVj  =  0 


:.4")  -^(Po-ro)a,p-  ^(po,Vo)ajafcVi  =  o. 

dp  d\ ki 

Consider  e.xponential  solutions  of  the  form 


^  =  ejp(i(x,0  +  At)  ^ 


where  (r.O  =  Substituting  2.5  into  2.3*-2.4*,  we  derive  a  generalized  eigenvalue 

problem  for  A: 


V  0 


Here  the  principal  part  of  5(0  block  structure 


0 

■f'sO  Q(0 


where  denotes  the  transpose  of  tho  column  vector  C  Performing  a  similarity  transfor¬ 
mation  to  eliminate  complex  entries.  /l  O  takes  the  form 


f.(s")  =  P0(V^)sS 


dp 
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and  Q  is  given  by 

(2.9) 

We  cull  the  Equations  2. 1-2.2  (linearly)  stable  if  the  real  part  of  A  is  negative:  Re  A(^)  <  0. 

Our  first  results  are  in  two  space  dimensions  and  relate  the  growth  of  the  generalized 
eigenvalue  A  to  the  eUipticity  of  the  steady  state  equations  associated  with  2.3*-2.4* 

(2.10)  PodjVj=0 

(2.11)  =  0 

Call  2.10-2.11  elliptic  if  detS{^)  ^  0  for  all  nonzero  ^  G  R^.  Roughly  speaJdng,  the 
theorem  below  states  that  Equations  2. 3-2. 4  become  unstable  when  2.10-2.11  change  type 
from  elliptic  to  hyperbolic. 

Theorem  2.1.  In  two  dimensions.  Equations  2. 3-2.4  are  stable  if  and  only  if  (i)  Equations 
2.10-2.11  axe  elliptic:  and  (ii)  both  eigenvalues  of  the  stress  tensor  T  are  non-negative. 

(Remark:  Shearer  and  Schaeffer  [14]  have  shown  that  the  results  of  this  theorem  remain 
true  for  the  full  system  2. 1-2.2.) 

One  implication  of  this  theorem  is  that  Equations  2. 3-2. 4  may  be  unstable  even  when 
the  stresses  are  on  the  consolidation  side  of  the  yield  locus  (i.e.  stress  states  for  w^hich 
tr{V)  >  0)  [9].  This  statement  contradicts  a  widely  held  view  that  deformation  is  stable  on 
the  consolidation  side  of  the  yield  locus  and  unstable  on  the  expansion  side  (i.e.  tr{V)  <  0); 
the  results  of  [9]  show  that  the  Critical  State  (i.e.  tr{V)  =  0)  is  the  boimdary  between 
well  pose  ]  and  ill  posed  behavior.  A  second  consequence  of  the  theorem  is  that  all  Fourier 
modes  along  the  specific  unstable  direction  lose  stability  simultaneously.  Which  specific 
mode  eventually  dominates  grovcth  (in  the  linear  regime)  depends  on  initial  and  boimdarj' 
conditions;  this  suggests  extreme  sensitivity  of  experiments  to  imperfections  in  laboratory 
protocol. 

The  analogue  of  Theorem  2.1  holds  in  three  dimensions,  but  the  geometric  interpretation 
is  more  complicated  [10].  In  particular,  when  the  steady  state  equations  lose  eUipticity 
they  do  not  become  hj-perbolic:  rather  they  are  of  mixed  type.  Because  of  the  complicated 
geometry  of  three  dimensional  flow.  .M.ACSVMA  is  being  employed  to  perform  the  algebra 
of  the  problem. 

§3.  NUMERICAL  SIMULATIONS 

Preliminary  efforts  at  numericaliy  integrating  Equations  2. 1-2.2  were  reported  in  the 
original  proposal.  Here  we  describe  our  current  work  on  this  phase  of  the  project. 

Consider  the  following  scenario  suggested  by  Theorem  2.1  for  the  development  of  shear 
bands.  Beginning  with  a  homogeneous  stable  flow,  material  consolidates  upon  deformation, 
eventually  reacliing  the  instability  locus.  Subsequent  deformation  remains  smooth,  but  it 
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is  no  longer  homogeneous  and  local  imperfections  grow  in  time.  At  some  point  in  the 
sample,  an  imperfection  ultimately  becomes  ill  posed  as  it  reaches  the  Critical  State  locus. 
Such  a  point  could  act  as  the  seed  of  a  velocity  discontinuity  which  then  expands  into  a 
shear  band.  Numerical  simulations  should  be  able  to  confirm  this  conjectured  route  to 
shear  band  formation,  but  the  ill  posedness  of  two  dimensional  deformations  precludes  the 
simple  intergration  of  Equations  2. 1-2.2.  We  follow  a  suggestion  of  Jackson  and  regularize 
the  system. 

Let  us  elaborate.  Near  a  shear  band,  there  is  large  momentum  transfer  over  distances 
on  the  order  of  a  few  particle  diameters.  Jenkins  [4],  Savage  [5],  and  colleagues  have 
developed  a  theory  for  granular  flow  at  very  high  sheair  rates  based  on  an  extension  of 
standard  kinetic  theory  to  slightly  inelastic  particles.  Incorporating  the  Enskog  dense  gas 
correction,  the  granular  continuum  equations  are  essentially  the  Navier-Stokes  equations 
for  a  compressible,  heat  conducting  fluid.  The  “heat'’  is  not  molecular  energy,  but  the 
so-called  “granular  temperature”  mecisuring  a  particle’s  deviation  from  the  mean  motion 
of  its  neighbors.  There  is  one  non-classical  term  in  the  theory,  a  heat  sink  which  represents 
the  energy  loss  due  to  inelastic  collisions.  Let  the  stress  tensor  derived  from  this  kinetic 
theory  be  T*,  and  denote  the  stress  tensor  of  the  Critical  State  theory  as  .  Jackson  [6] 
considers  a  system  consisting  of  Equations  1.1-1. 2  and  the  energy  equation  of  the  kinetic 
theory,  with  the  total  stress  being  the  sum 

(3.1)  T  =  T^  +  T^-, 

^6:  also  shows  how  to  derive  appropriate  boundary  conditions  for  flow. 

The  primary  consequence  of  the  addition  of  these  kinetic  stresses  is  to  add  an  overall 
regularization  to  the  system  2. 1-2.2.  In  unpublished  notes,  Schaeffer  and  Pitman  demon¬ 
strate  that  the  system  1.1-1. 2.  3.1  is  well  posed.  However,  near  unstable  and  iU  posed 
states  the  system  possesses  only  a  small  amount  of  damping  of  very  high  frequency  Fourier 
components.  For  moderate  wave  numbers  near  these  states  there  is  a  competition  between 
the  growth  of  the  friction  terms  and  the  damping  of  the  kinetic  terms. 

In  order  to  numerically  solve  the  equations  without  extreme  time-step  limitations  and 
to  adequately  account  for  the  potential  growth  of  Fourier  modes  as  outlined  above,  we  use 
an  implicit  Beam- Warming  type  scheme  [1]  to  integrate  the  equations.  Beam- Warming 
allows  for  extension  to  fully  multi- dimensional  computations.  The  hyperboUc  terms  are 
treated  seperately  from  the  parabolic-like  terms.  Following  Harten  and  Yee  [15]  the  hy¬ 
perbolic  terms  are  differenced  in  -'n-h  a  way  as  to  (try  to)  preserve  the  Total  ^’ariation 
Diminishing  character  of  such  tt  r;;.'.  Parabolic-like  terms  are  centered  differenced,  and 
the  entire  scheme  is  constructed  m  d'  ira"  form.  The  resulting  set  of  difference  equations 
involve  block-tridiagonal  matri.x  : '•.  "u.  Such  schemes  have  been  successfully  applied  to 
aerodynamic  calculations. 

Our  first  numerical  experim'  -  •  •  ;. fined  to  quasi-one-dimensional  flow  in  a  hopi^er. 

See  the  attached  figure,  where  v.  -  .•  •  the  velocity.  a.s  a  fiinction  of  position,  at  various 

times.  The  geometry  of  this  {>:  ■  -imilar  to  that  for  a  spherically  expanding  wave 

passing  through  the  ground,  and  formulation  is  simpler  in  the  hopper  problem. 

In  the  simulation,  we  consider  inua-.dy  ;,dy  flow;  at  time  /  =  0.  we  suddenly  increased 
the  flow  rate  at  the  exit.  We  have  examin'd  the  transient  behavior  as  flow  sc'ttles  into  a 
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new  quiescent  steady-state  regime.  Velocity  waves  propagate  quickly  through  the  material 
and  equilibrate  rapidly;  the  solid  line  in  the  figure  is  the  steady  state,  attained  after  about 
10  seconds.  Note  the  overshoot  of  the  velocity  profiles.  Density  vairiations  remain  active 
for  an  extended  period  of  time  and  Stresses  become  negative  during  flow,  signaling  tensile 
loads.  In  many  granular  flow  simulations,  tensile  stresses  lead  to  catcistrophic  numerical 
instabilities;  the  visco-plastic  formulation  (3.1)  allows  us  to  continue  computations  through 
the  tensile  region. 
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